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L Introduction 


(U) A combined Monte Carlo — CFD algorithm has been developed re- 
cently at NASA lewis Research Center for turbulent reacting flows. In this al- 
gorithm, conventional CFD schemes are employed to obtain the velocity field 
and other velocity related turbulent quantities, and a Monte Carlo scheme 
is used to solve the evolution equation for the probability density function 
(pdf) of species mass fraction and temperature. 

(U) In combustion computations, the predictions of chemical reaction 
rates (the source terms in the species conservation equations) are poor if 
conventional turbulence models are used. The main difficulty lies in the fact 
that the reaction rate is highly non-linear, and the use of averaged tempera- 
ture produces excessively large errors. Moment closure models for the source 
terms have attained only limited success. The probability density function 
(pdf) method seems to be the only alternative at the present time that uses 
local instantaneous values of the temperature, density, etc., in predicting 
chemical reaction rates, and thus may be the only viable approach for more 
accurate turbulent combustion calculations. The closure problem and the 
need for pdf has been discussed in detail in a previous paper[l]. 

(U) Assumed pdf are useful in simple problems; however, for more gen- 
eral combustion problems, the solution of a evolution equation for the pdf is 
necessary. 


II. Method 

(U) Conventional CFD flow solvers axe quite common and needed no 
explanation here. We will concentrate on the derivation and solution of the 
pdf evolution equation in this section. 

2.1 PDF evolution equation of mass fraction. 

(U) For simplicity, the pdf equation for a single scalar will be derived 
here. The extension to multi species is trivial. The species transport equa- 
tions ( for a single species ) can be written as 

pd t Y = -pu k dkY + pdk(DdkY) + pw (0.1) 

where it* is the velocity, Y is the mass fraction, and w is the chemical source 
term; summation for repeated indices is understood. 
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Let P(Y) be the probability density function, and define <j> = eip(:Ac), 
then the characteristic function of the pdf can be written as the ensemble 
average of <j>: 

< 4> >= / P(Y)4>dY, 

and the pdf can be written as the inverse Fourier transform of the character- 
istic function: 

P(Y) = ^j<<t>> exp{-iXY)dY. 

Now differentiate <f> with respect to time, 

d t <j> = i\exp{i\Y)(d t Y), 

and replace dtY in the above equation with the right hand side of eq. (1), 
we obtain 

pd t <i> = i\exp(i\Y)[-pu k d k Y + d k (pDd k Y ) + pw], 

which can be rearranged as 

d t (p<j>) -f d k (pu k 4>) = i\<f>d k (pDd k Y) + i\<f>pw 

Take ensemble average of the above equation gives 

dt< p<t>> +d k < pu k <f> >= i\ < pw<f> > +*A < <f>d k (pDd k Y) > . 

The inverse Fourier transform of the above equation gives the evolution equa- 
tion for the probability density function: 

d t (pP) -1- d k (p <u k > P) + dy(pwP) 

= -d k (p < u' k \ Y >P)~ dy{< d k ( P Dd k Y)\Y > P) 

where the terms represent mean convection, chemical reactions, turbulent 
convection, and molecular mixing, respectively. 

For multi species and temperature, the pdf is written as P(ip lt ..., V'n), 
where represent the scalars including mass fractions and tem- 

perature. A similar derivation gives 

N 

pdtP + pv a d a P + p'%2dj i {w i (iJ> ly ...,ilf r f)P} 

t=l 

= ~d a (p < t£| rf>i > P) - > -P)- 

i=i ,=i 
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(U) The left hand aide of the above equation can be evaluated exactly 
and requires no modeling; the right hand side terms contain the conditional 
average of the Favre velocity fluctuation and the conditional average of the 
scalar dissipation and require modeling. 

2.2 Modeled pdf equation 

(U) The first term that needs modeling is the turbulent diffusion term, 
for which the following gradient model is used: 

- < u ' k |&- > P & D t d a P , 

where D t is the turbulent diffusion coefficient, which is set to be equal to the 
eddy viscosity, i.e., the turbulent Schmidt number is equal to one. 

(U) The next term needs modeling is the molecular mixing term. A 
coalescence/dispersion model is used for this term, which has the following 
general form: 

-pYl'52 d l i 4j( < > P ) 

•«i i=i 

-T'J / pw)p(r)Tm\rwdj>" - p 

= D m P 

where T is the transition probability. For more detail description of the 
molecular mixing model, see Ref 2. The modeled the pdf equation is then 

N 

pd(P + pv a d a p + p^2 d^{wi{rl> 1 , ..., rJ)ff)P} 

»=1 

= d j (D t d a P) + D m P 


2.3 Solution of the pdf equation 

(U) A fractional step method is used to solve the pdf equation in a 
Monte Carlo simulation. Three processes, namely, (1) convection/diffusion, 
(2) molecular mixing, and (3) chemical reaction, are simulated consecutively: 

N 

+ - D m )(d t + pvd a -d a (D t d a )P = 0. 

<=i 


71 


<u V ’>/U 


This equation is solved using a Monte Carlo simulation, i.e., the pdf is 
represented by an ensemble of samples. These sample points move in the 
hyperspace of (xi,X2,X3,^>i,...,^ n ) according to the above equation. 


III. Results 

(U) A heated turbulent plane jet had been calculated and results re- 
ported at last NASP symposium. The case has since been recalculated using 
modified turbulence model, and improved results are shown here. 


Velocity distribution 




Turbulent shear stress 
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(U) A hydrogen— fluorine diffusion flame was calculated. The flowfield 
and flow conditions, as well as the calculated velocity and turbulent kinetic 
energy distribution, are given on this page. 



Velocity distribution 


Turbulent kinetic energy 




(U) Figures show the calculated heat release of the H2-F2 flame (tem- 
perature distribution) compared to experimental data, and the species mass 
traction from pdf calculation. 






(U) Computed pdf of temperature distribution in the H2-F2 flame. 
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